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ABSTRACT 

We  investigate  the  small  sample  behavior  and  convergence  properties  of  confidence  interval 
estimators  (CIE’s)  for  the  mean  of  a  stationary  discrete  process.  We  consider  CIE’s  arising 
from  nonoverlapping  batch  means,  overlapping  batch  means,  and  standardized  time  series, 
all  of  which  are  commonly  used  in  discrete-event  simulation.  For  a  specific  CIE,  the 
performance  measures  of  interest  include  the  coverage  probability,  and  the  expected  value 
and  variance  of  the  half-length.  We  use  both  empirical  and  analytical  methods  to  make 
detailed  comparisons  regarding  the  behavior  of  the  CIE’s  for  a  variety  of  stochastic 
processes.  All  of  the  CIE’s  under  study  are  asymptotically  valid;  however,  they  are  usually 
invalid  for  small  sample  sizes.  We  find  that  for  small  samples,  the  bias  of  the  variance 
parameter  estimator  figures  significantly  in  CIE  coverage  performance  -  the  less  bias  the 
better.  A  secondary  role  is  played  by  the  marginal  distribution  of  the  stationary  process.  We 
also  point  out  that  not  all  CIE’s  are  equal  -  some  require  fewer  observations  before 
manifesting  the  properties  for  CIE  validity. 


Subject  classifications:  Simulation,  statistical  analysis,  statistical  estimation,  time  series, 
small  sample  behavior  of  confidence  intervals 
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This  paper  studies  the  small  sample  behavior  and  convergence  properties  of 
a  number  of  confidence  interval  estimators  (CIE's)  for  the  mean  p  of  a 
stat ionary  process,  X^, — ,X^.  These  CIE's  are  typically  of  the  form 

PrJM  E  \  t  =  I  -  (1) 

where  X  =  Z.X./n,  t  .  ^  is  the  >-quantile  of  the  t-distribution  with  d  degrees 
n  1  i'  d.Y  ^  *• 

of  freedom  (d.o.f.)f  and  V  estimates  =  nVar(X^)  (or  the  variance  parameter 

2  2  2  2 
a  =  lim  o  ).  A  "good"  estimator  for  o  (or  a  )  is  the  cornerstone  of  a 
n-»ra  n  ^  n^ 

valid  CIE  for  p.  Many  estimators  have  been  studied  in  the  context  of  discrete- 
event  simulation:  nonoverlapping  batch  means  (NOBU)  (Conway  1963,  Schmeiser 
1982,  Kang  1984);  independent  replications;  overlapping  batch  means  (OBM) 
(Meketon  1980,  Meketon  and  Schmeiser  1984);  standardized  time  series  (STS) 
(Schruben  1983,  Goldsman  1984,  Glynn  and  Iglehart  1990);  spectral  theory 
(Fishman  1973,1978,  Heidelberger  and  Welch  1981,1983);  ARMA  time  series 
modeling  (Fishman  1973,1978,  Schriber  and  Andrews  1984);  and  regeneration 
(Crane  and  Iglehart  1975,  Crane  and  Lemoine  1977,  Fishman  1978). 

There  is  considerable  work  which  compares  the  various  CIE  methodologies. 
The  Monte  Carlo  (MC)  work  mainly  deals  with  small  sample  CIE  performance;  see, 
e.g..  Law  and  Kelton  (1984)  and  Goldsman,  Kang,  and  Sargent  (1986).  Analytical 
results  are  almost  all  asymptotic:  Goldsman  and  Schruben  (1984),  Goldsman  and 
Meketon  (1986),  Damerd.ji  (1987),  Glynn  and  Iglehart  (1990),  and  Schmeiser  and 
Song  (1989)  all  compare  various  combinations  of  the  CIE's. 

In  this  paper  we  study  finite  sample  behavior  of  CIE's  from  NOBM,  OBM,  and 
STS.  Section  1  gives  background  material.  Section  2  reports  on  statistical 
properties  of  various  variance  estimators.  Section  3  presents  analytical 
results  for  some  sprrial  cases  and  then  summarizes  a  MC  study  of  the  CIE's.  In 
Part  I  of  our  MC  work,  we  break  the  n  observations  into  b  batches,  and  we 
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observe  what  happens  as  the  batch  size  grows.  For  ''small"  b,  NOBM  performs  the 
best  with  respect  to  CIE  coverage.  For  "large"  b,  both  NOBM  and  OBM  fare  the 
best.  However,  an  STS  method  produces  CIE's  with  smaller  expected  lengths. 
Another  comparison  is  carried  out  in  Part  II  of  our  MC  work,  where  we  fix  the 
d.o.f.  d  and  oDserve  wliat  happens  as  n  grows;  here,  the  NOBM,  OBM,  and  STS 
"combined"  CIE's  perform  similarly.  Section  4  discusses  our  findings,  and 
Section  5  summarizes,  ffe  conclude  that  the  bias  of  V  is  the  most  important 
factor  in  determining  a  CIE's  validity;  a  secondary  role  is  played  by  the 
marginal  distribution  of  the  X^'s.  We  also  find  that  a  CIE  having  superior 
large  sample  properties  may  have  relatively  poor  small  sample  performance.  We 
offer  practical  and  research  recommendations. 


1.  BACKGROUND 

We  review  the  CIE's  and  stochastic  processes  under  study.  We  assume  the 

stochastic  processes  satisfy  certain  moment  and  mixing  conditions,  as  described 

2  2 

in  the  cited  papers.  We  henceforth  use  the  notation  Nor(p,T  ),  X  (d),  X(d), 
Exp(x),  and  t(d)  to  represent  the  normal,  chi-square,  chi,  exponential,  and  t 
distributions,  each  with  the  appropriate  parameters. 


1 . 1  Nonoverlapping  Batch  Means 

Suppose  we  divide  X^,...,X^  into  b  >  1  adjacent,  nonoverlapping  batches  of 

size  m  (assume  n  =  bm) .  The  i-th  batch  mean,  i  =  l,...,b,  is  X.  = 

'  ^  i,m 

y™  4  X/ .  ,  ./m.  In  implementing  the  method  of  NOBM,  we  assume  the  X.  's  are 

2  2 
approximately  i.i.d.  Nor(p,a  /m).  The  NOBM  estimator  for  o  is 

=  m  X.^^  -  ]^/(b-l)  $  a^X^(b-l)/(b-l), 


where  ”  -♦  "  denotes  convergence  in  distribution  as  m  -*  ®.  The  NOBM  CIE  for  p 
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is  given  by  (1)  with  d  =  b-1  and  V  = 

1.2  Overlappine  Batch  Keans 

Define  the  i-th  overlapping  batch  mean,  i  =  l,  —  ,n-in+l,  by  X(i,m)  = 
y™  n  X.,  ./m.  The  OBM  estimator  for  is 

^j=0  i+j' 

Vq  =  nm  [X(i,m)  -  l^]^/[(n-m+l)(r-m)]  , 

and  is  almost  identical  to  Bartlett's  spectral  estimator  (see  Priestley  1982). 

The  OBM  CIE  for  p  is  given  by  (1)  with  V  =  V^;  its  validity  depends  on  being 
2  2 

approximately  a  X  (d)/d.  Meketon  and  Schmeiser  (1984)  take  d  =  1.5* (b-1),  where 
b  =  n/m.  Based  on  MC  experimentation,  Schmeiser  (1986)  recommends  d  = 
1.5*(b-l)[l+(b-l)  we  shall  use  this  value  in  our  MC  work. 


1.3  Standardized  Time  Series 

Suppose  Xj, — ,X^  is  divided  into  b  ^  1  adjacent,  nonoverlapping  batches 

of  size  m.  For  i  =  l,...,b,  let  A.  =  [(m+l)/2  -  .  Schruben 

(1983)  assumes  the  A^'s  are  approximately  i.i.d.  normal,  and  proposes  the  area 

2 

and  combined  estimators  for  a  ; 


V  =  A^ 

A  /  3  ^1=1  1 

(m  -m)b 


2v2/.  ^ 

o  X  (bj 


(b-l)V,  +  bV^ 


o^X^(2b-i; 


CIE's  for  p  are  formed  by  substituting  the  appropriate  V  and  d  in  (1).  Schruben 

2 

also  derives  the  so-called  "maximum"  estimator  for  a  (see  Subsection  4.2). 


1.4  Some  Time  Series  Processes  of  Interest 

We  will  have  occasion  to  use  the  following  ARMA-type  processes. 


MA(l):  X^  =  where  ~  i.i.d.  Nor(0,l)  and  -1  <  9  <  1. 

MA'(l):  X.  =  e.  +  9e.  ,,  where  e.  ~  i.i.d.  Exp(l)  and  -1  <  9  <  1. 

AR(1):  X^  =  where  £.  ~  i.i.d.  Nor  (0,  l-tp^) ,  -l<tp<l. 


-4- 


EAR(l):  X.  =  ^ 


<pX 


i-l 


w.p.  (p 


I  +  Cj  w.p.  l-(p 


where  ~  i.i.d.  Exp(l)  and 
0  ^  ip  <  1  (see  Lewis  1980). 


We  also  consider  other  stationary  normal  processes  as  well  as  the  waiting  time 


(delay)  processes  of  U/M/1  and  E^/K/l  queues. 


1.5  CIE  Performance  Criteria 

Denote  the  NOBM,  OBM,  area,  and  combined  CIE's  by  CIEj^,  CIE^,  CIE^,  and 
CIEj,,  resp.  The  half-length  of  a  generic  CIE  is  H  =  t^ 

the  following  CIE  performance  measures:  coverage  (Pr|p  £  i  H0>  B[HJ.  and 
Var(H).  Among  CIE's  which  achieve  coverage  1  -  a,  we  prefer  that  with  the 
smallest  E[H],  and  then  that  with  the  smallest  Var(H). 


2.  PROPERTIES  OF  VARIOUS  VARIANCE  ESTIMATORS 

We  give  some  results  on  the  bias  and  variance  of  the  NOPy,  OBk,  and  STS 
estimators,  and  on  the  asymptotic  performance  of  the  corresponding  CIE's. 

2 . 1  Bias  of  the  Estimators 

The  bias  of  V  as  an  estimator  for  o  is  Bias(V)  =  a  -  E[Vj.  Goldsman  and 
keketon  (1986)  show  that  lim  m lim,  Bias(V„)  =  lim  m  lim,  Bias(V_)  = 

ID-tco  i>+®  N  IIh>®  l>-»®  0 

lim^^m  limj^^Bias(V^)/2  =  1  im^^m  Bias(V^)/3.  All  of  these  estimators  are 
asymptotically  unbiased  as  m  ®,  but  for  small  m,  these  estimators  can  be 
quite  biased. 

2 

Example  1:  For  the  AR(l)  and  EAR(l),  the  Appendix  gives  a  =  (l+ip) /(l-ip)  and 

E[VJ  =  ,  ^2  .  2v 

mb(l-(p)^  m(b-l)  (l-(p)^  m(l-(p)^ 


for  large  m  and  b,  (2) 
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E[V.]  =  a  - 


m(b-l)(l-(p) 


.  mb  T  ,  1 1  / «  mo-mT 

_  1  - ^  I  +  _ 

l-(p)^  b^  J  m(b-l)  (mb-m+1)  (1- 


lUv/.  mb-m+ls 


(b-1)  (mb-m+1)  (l-(p)' 


.  2 

=  a  - 


m(l-(p) 


■g  for  large  m  and  b. 


E[V  ]  =  +  3  ^ -"f - g 


^+l-(p°*(m+l)^  .  (p-(m+l-mq 


r  ^ ,  4  I 
-m  +l~(p 
<  » 

4 


2  6u 
a  -  "  J 


m(l-ip) 


for  large  m. 


E[V^]  =  0*"  + 


(2b-l)m(l-(p) 


,  mb 

4b  +  i:^ 


.  2 

=  a  - 


m(l-(p)‘ 


for  large  m  and  b. 


If  m  and  b  are  large,  the  bias  results  anticipated  by  Goldsman  and  lleketon  are 
attained.  Table  1  contains  exact  E[V]'s  for  b  =  2  and  16  and  various  m.  For 
small  b,  Bias(Vj^)  <  BiasCV^)  <  Bias(V^)  <  BiasCV^);  for  large  b,  Bias(Vjj)  = 
BiasCV^)  <  Bias(V^)  <  Bias(V^).  The  biases  become  negligible  as  m  grows. 

Example  2:  For  the  MA(l)  and  MA'(l),  the  Appendix  gives  =  (1+0)^  and 

E[Vjj]  =  -  20(b+l)/mb  =  -  2H/m  for  large  b.  (6j 


20  rb  +1 


■  ^  ,  ,  I  +  -  20/m  for  large  b. 

mb-m+1 


E[V^]  =  a  -  60/m. 

E[V^]  =  -  (4b+2)0/mb  =  o^  -  40/m  for  large  b. 

The  conclusions  from  Example  1  again  hold.  ,, 


Although  Bias(V)  is  interesting  in  its  own  right,  the  bias  directly 
affects  CIE  performance.  Consider  the  unrealistic  case  in  which  m  is  fixed  and 
b  -»  CD.  For  the  estimators  studied  here,  one  can  show  that  as  b  -»  <d,  V  -  E[V] 
w.p.l,  and  so  T  =  (l^-p)/(V/n)'^  2  Nor (0,a^/E[V]) .  Falsely  assuming  T  ~ 

Nor(0,l)  as  b  00  yields  incorrect  CIE  coverage. 
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[zi_^/2(E[V]/a2)^]  -  1  .  (10) 

where  ♦(•)  and  are  the  Nor(0,l)  c.d.f.  and  Y-quantile.  If  E[V]/a^  <  [>]  1, 
then  coverage  <  [>]  1  -  a.  (Coverage  is  quite  sensitive  to  decreases  in 

E[V]/a  .)  So  the  less  bias  the  better.  As  b  -*  ®,  CIE„  and  CIE^  tend  to 

achieve  the  desired  coverage  more  quickly  with  respect  to  m  than  <^o  CIE^  and 
CIE^;  see  Sargent,  Kang,  and  Goldsman  (198/)  (S-K-G) . 

2.2  Variance  of  the  Estimators 

Goldsman  and  Meketon  report  that  as  m  and  b  become  large,  b*Var(Vj|^)  -*  2a  , 
b*Var(VQ)  4a‘^/3,  b-Var(V^)  -»  20“^,  .nd  b-Var(Vj,)  (cf.  Damerdji  1987). 

For  i.i.d.  Xj,...,X^,  Kang  and  Goldsman  (1990)  find  Var(Vj^)  and  Var(V^) 
exactly,  as  do  Song  and  Schmeiser  (1989)  for  Var(VQ).  Exact  results  for  other 
processes  and  for  Var(V|^)  are  tedious  to  derive.  Of  course,  one  can  also 

calculate  the  mean  squared  errors  of  the  V’s  (cf.  Schmeiser  and  Seng  1989). 


2.3  Asymptotic  Properties  of  the  CIE's 

Schmeiser  (1982)  and  Goldsman  and  Schruben  (1984)  note  that  as  m  -♦  », 

(mb)^'H  5  <rt^^j_„/2X(d)/>ld  (11) 

for  CIE||^,  CIE^,  and  tJIE^.  An  analogous  approximate  result  holds  for  CIE^^. 

Under  (11),  the  CIE's  achieve  coverage  1  -  «  as  m  -♦  ®.  Further,  if  r(*)  is  the 
gamma  function,  then 

/  Till  /'.^/■\Xr((d+l)/2)  , 


mbVar(H)  -♦  ^^t^ 


.l-a/2  ( 


r(d/2) 

and 

(12) 

r((d+l)/2)  ' 
r(d/2) 

D- 

(13) 

d  and,  hence. 

in  b. 

So  for  large 

m  with  fixed  b,  E[Hj^]  >  E[H^]  >  E[HqJ  >  ELH^.]  and  Var(Hj^)  >  Var(H^)  >  Var(HQ)  ^ 


Var(H^),  the  subscripts  having  the  obvious  meanings.  Goldsman  and  Schruben 
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(1984)  and  Meketon  and  Schmeiser  (1984)  b  ®  in  (12)  and  (13)  to  get 


lim 


E[H] 


E[H.] 


1  >  H  -  Hq,  , 


and 


lim 


b-»® 

n»4ffl 


Var(H) 


Var(Hj^) 


'  1.  H  =  H^ 

.  2/3,  H  =  H^j  . 

.1/2.  H  =  H^ 


(14) 


Thus,  as  b  also  becomes  large,  all  of  the  CIE’s  have  about  the  same  E[H]'s,  but 
CIE^  has  the  smallest  Var(H). 


3.  rlNITE  SAMPLE  CONFIDENCE  INTERVAL  ESTIVIATIC.. 

Small  sample  analysis  of  CIE's  is  difficult-  We  present  a  few  exact 
results,  but  most  of  the  section  is  devoted  to  a  UC  study. 

3. 1  Some  Analytical  Results 

Example  3:  Suppose  X^,...,X^  --  i.i.d.  Nor(|i,^').  Then  ~  T^X^(b-l)/(b-l) , 

V,  ~  T^X^(b)/b,  and  V„  ~  T^X^(2b“l)/(2b-l) .  Further,  X  is  independent  of  V^, 

A  L  n  iw 

V^,  and  (see  Appendix).  So  (l)  holds  exactly  for  CIEj^,  CIE^,  and  CIE^,.  We 
could  not  obtain  such  results  for  CIE^  or  for  nonnormal  i.i.d.  processes, 

Example  4:  We  can  deiive  exact  results  for  CIE.  when  b  =  1  (n  =  m)  and  X  ,...,X 

A  J.  O 

is  stationary  normal.  Then  ~  Nor  (0,E’ A^]) ,  V  ~  E[V  ]X^(l),  and  X  is 

1  1  A  A  n 

normal  -^nd  independent  ^f  V  (see  Appendix).  So  (X  -u)(nc/V  )'^  ~  t(l),  where 

A  HA 

c  E  Hence,  the  coverage  of  CIE^  is  2Pr5t(l)^tj  ” 

t2/K)Tan  ^(t^  l-a/2'^*^^'  (10)»  the  coverage  is  sensitive  to  decreases 

in  c.)  Similarly.  E[H.l  -  t.  ,  /o(2E[V, 1/nn)^.  To  illustrate,  suppose  the 

A  1 ,  1~0(/ d  A 

X.'s  are  AR(l).  Figure  1(a)  uses  Example  1  and  (A-3)  to  plot  coverage  vs. 
log^n  for  1  -  O'  =  0.90  and  various  ip.  We  see  that  for  ip  >  [<]  0,  the  coverage 
is  <  [>]  1  -  a.  As  |(p|  approaches  0  or  as  n  grows,  Bias(V^)  decreases,  and  the 
coverage  approaches  1  -  or.  Figure  1(b)  has  analogous  plots  of  E[H^J  vs. 
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2  “X 

loggii.  If  (p  =  0.0,  then  E[V^]  =  a  =1,  and  so  E[H^]  decreases  at  rate  n 

If  (p  =  -0.9,  then  E[V^]  decreases  to  a  =  1/19,  and  E[H^]  decreases  at  rate  n 

(after  initially  decreasing  faster).  The  (p  =  0.9  plot  for  E[H^]  increases  and 

then  decreases  as  n  grows.  This  occurs  since  E[V^]  increases  to  o  =19  as 

n  ®,  while  the  -Jn  in  E[H^]'s  denominator  becomes  large,  jj 

Example  5:  We  give  exact  results  for  CIBj^  when  b  =  2  (n  =  2m)  and  Xj^,...,X^  is 

stationary  normal.  Then  V  =  m2.(X.  -1  )^/(b~l)  =  m(X,  -  X„  )^/2;  so  V„  ~ 

N  1  1  ,m  n  1  ,m  2,m  N 

""  2  = 

E[Vj^]x  (l).  Since  X^  is  normal  and  independent  of  for  b  =  2  (see  Appendix), 
we  have  (X^-p) (nc ’/Vj^)^  ~  t(l),  where  c'  =  E[Vj^]/a^.  The  coverage  is 
(2/jt)Tan  ^(t^  ^i  ® 

AR(1),  then  (A-3)  yields  c^,  and  (2)  with  b  =  2  gives  E[Vj^];  c',  coverage,  and 
E[H|^]  then  follow.  These  performance  measures  behave  as  in  Example  4. 

It  is  difficult  to  generalize  the  above  CIE  results  to  b  >  2,  since  we 

2 

would  then  have  correlated  X  random  variables.  Thus,  we  only  gave  exact 
results  for  simple  cases.  We  resort  to  MC  experimentation  in  the  sequel. 

3.2  Desi£n  of  the  Monte  Carlo  Study 

Our  goal  was  to  assess  CIE  performance  over  a  variety  of  stochastic 
processes  and  choices  of  number  of  observations  n,  batch  size  m,  d.o.f.  d,  and 
desired  coverage  1  -  a.  We  simulated  the  following  processes:  AR(l)  with 
(p  =  0.0, iO.l, +0.5,10.9;  EAR(l)  with  tp  =  0.0. 0.1, 0.5, 0.9;  MA(l)  and  MA'(l)  with 
0  =  +0.1,10.9;  M/M/1  with  traffic  intensity  p  =  0.6, 0.8  (service  rate  =  1.0); 
and  Eg/M/1  with  p  =  0.6  (service  rate  =  1.0).  Each  run  was  initialized  from  the 
appropriate  steady  state  distribution.  All  uniform  [normal]  variates  were 
generated  from  algorithm  UNIF  [TRPNRM]  in  Bratley,  Fox,  and  Schrage  (1987); 
exponentials  used  inversion. 
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For  ease  of  exposition,  we  divided  the  study  into  two  parts.  In  Part  I  we 

fixed  1  -  a  and  b,  and  then  charted  CIE  coverage  as  a  function  of  m.  Roughly 

speaking,  we  wanted  to  know  which  CIE  first  achieved  acceptable  coverage  as  m 

increased  with  fixed  b.  Further,  which  CIE’s  had  the  best  E[H]  and  (to  a 

lesser  extent)  Var(H)?  For  the  stochastic  processes  discussed  above,  we 

conducted  at  least  1000  independent  runs  of  16384  observations;  these  runs  were 

used  to  calculate  CIE.,,  CIE^,  CIE.,  and  CIE„  and  the  resulting  performance 

W  U  A 

characteristics  for  all  choices  of  m  =  2  ,  k  =  0,' . 10,  b  =  1,2,4,8,16,  and 

1  -  a  =  0.80,0.90,0.95,0.99. 

In  Part  II,  we  set  1  -  a  =  0.90  and  d  =  3  and  15,  and  then  charted 
coverage  as  a  function  of  n.  For  d  =  3  [d  =  15],  CIEj^  uses  b  =  4  [b  =  16], 

CIEq  uses  b  =  2  [b  =  11],  CIE^  uses  b  =  3  [b  =  15],  and  CIE^,  uses  b  =  2 
[b  =  8].  We  conducted  our  experiments  on  a  number  of  stochastic  processes, 
each  of  which  used  2000  runs  of  16384  observations  to  calculate  the  CIE's  for 
all  n  =  2*^,  k  =  4,..., 14,  and  d  =  3  and  15.  Since  d  is  fixed,  (ll)-(13) 
suggest  that  all  of  the  CIE's  will  perform  about  the  same  for  large  m;  however, 
we  suspected  that  they  would  exhibit  different  small  sample  performance  because 
the  variance  estimators  incorporated  in  the  CIE's  operate  under  different 
assumptions.  We  first  discuss  the  underlying  assumptions  for  Vj^,  V^,  and 
since  these  estimators  require  independence  between  batches  (V^  does  not).  The 
estimator  V„  [V.]  assumes  that  the  X.  's  [A.'s]  are  i.i.d.  normal;  the 
combined  estimato,  must  satisfy  both  assumptions.  For  fixed  d  and  n,  and 
use  roughly  half  the  batch  size  of  V^;  hence,  the  assumption  of  i.i.d. 

X  's  [A.'s]  is  harder  to  achieve  for  V„  [V.]  than  for  V„.  On  the  other  hand, 
V„'s  assumption  of  normality  of  the  X.  's  is  easier  to  satisfy  than  V  's 
assumption  of  normality  of  the  A.'s  which  relies  on  a  more  restrictive 
functional  central  limit  theorem.  The  normality  question  for  vs.  is  not 
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as  clear  since  V„  uses  half  of  V_'s  batch  size.  CIE_  appeals  to  spectral 
n  L  U 

^  2  2 

theory  co  directly  assume  that  is  a  X  (d)/d;  this  supposition  is  not  true 
for  nonnormal  processes  or  for  finite  batch  sizes. 

3.3  Results  from  Part  I  of  the  Monte  Carlo  Study 
3.3.1  Representative  results 

We  discuss  typical  results  from  Part  I.  Figures  2,  3,  and  4  illustrate 

CIE  performance  when  X^,  —  ,X^  are  AR(l)  ((p  =  0.9),  EAR(l)  ((p  =  0.9),  and  M/M/l 

(p  =  0.8),  resp.,  and  1  -  a  =  0.90.  The  AR(l)  and  EAR(l)  have  the  same 

covariance  functior  but  the  AR(l)  has  normal  marginals  while  the  EAR(l)'s  are 

exponential.  The  M/M/1 's  joint  distribution  is  more  complicated.  We  only 

consider  the  cases  b  =  2  and  16  since  CIEj^,  CIE^,  and  CIE^  require  b  ^  2,  and 

since  one  can  argue  that  b  =  16  is  "large"  (at  least  for  the  "usual"  choices 

of  a).  Each  of  Figures  2,  3,  and  4  has  four  graphs:  (a),(b)  are  for  b  =  2,  and 

(c),(d)  are  for  b  =  16.  In  (a)  and  (c),  we  plot  the  sample  coverage  (CVG)  vs. 

loggffl;  (b)  and  (d)  plot  the  sample  E[H]  (EHL)  vs.  log^m.  The  standard  error  of 

any  CVG  is  about  [CVG- (l-CVG)/1000]^. 

Figure  2  is  for  the  AR(1)  with  (p  =  0.9.  All  CVG's  are  poor  for  small  m, 

but  approach  1  -  a  =  0.90  as  m  increases.  For  b  =  2,  the  CVG  of  CIEj^  is 

closest  to  0.90;  for  b  =  16,  both  CIE„  and  CIE_  yield  the  best  CVG's.  This 

n  u 

makes  sense  since  V„  (for  b  =  2  and  16)  and  V„  (for  b  =  16)  are  less  biased 

N  u 

than  V^  and  V^^  (Example  1  and  Table  1).  A  related  consequence  is  that  CIEj^  and 
CIEq  produce  larger  EHL's  than  those  of  CIE^.  Note  that  the  EHL's  in  Figure  2 
increase  and  then  decrease  as  m  increases  -  the  same  bias-related  pattern  as  for 
the  exact  E[H^]  in  Figure  1.  The  EHL's  for  b  =  16  and  m  ^  128  from  Figure  2(d) 
are  more  or  less  in  agreement  with  the  limiting  E[H]/E[Hj^]  =  1  ratios  in  (14); 
this  is  one  reason  why  we  regard  b  =  16  as  "large." 
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Table  2  has  CVG's  for  the  AR(l)  with  b  =  16,  1  -  a  =  0.90,  and  various 
m  and  (p.  These  roughly  agree  with  the  b  eo  results  from  (10).  We  see  that 
CVG  <  [>]  0.90  when  ip  >  [<]  0.0.  For  any  (p  and  fixed  m.  Table  2  shows  that  the 
CVG's  of  CIE„  and  CIE„  are  closer  to  1  -  a  than  those  of  CIE,  and  CIE_. 

W  U  A  L 

Results  for  the  EAR(l)  process  with  (p  =  0.9  are  found  in  Figure  3,  whose 

plots  bear  a  striking  resemblance  to  those  of  the  AR(l)  in  Figure  2.  The  only 

notable  difference  between  the  two  figures  is  that,  for  fixed  m,  the  EAR(l)  has 

smaller  CVG's  (see  Table  1).  Since  the  AR(l)  and  EAR(l)  have  the  same 

covariance  structure,  the  EAR(l)'s  poorer  CVG's  are  probably  due  to  its 

exponential  marginals.  It  seems  that  the  effect  of  Bia8(V)  on  coverage  is  more 

significant  than  that  of  the  marginals  of  the  X^'s. 

Figure  4  concerns  the  M/M/1  process  with  p  =  0.8.  Again,  the  patterns  in 

the  figure  are  not  much  different  than  those  of  the  AR(1)  and  EAR(l).  The 

M/M/1  simply  requires  more  observations  to  attain  valid  coverage.  The  positive 

2 

serial  correlation  of  the  M/M/1  causes  the  estimators  for  a  to  be  biased  too 
low;  so  poor  coverage  results  for  small  m. 

3.3.2  Additional  Part  I  results 

In  S-K-G,  we  also  give  detailed  discussions  on  (among  other  things): 

-  CIE  performance  for  the  MA(l)  as  9  varies.  We  find  that  the  CIE's 
qualitatively  perform  about  the  same  as  those  for  the  AR(1). 

-  The  sample  •JVar(H)  (SHL)  performance  measure.  For  small  m,  the  SHL's  exhibit 
the  same  general  behavior  as  the  EHL's;  as  m  and  b  become  large,  the  SHL's 
behave  as  in  (14).  Even  though  the  ratios  from  (14)  are  manifested,  the 
differences  between  SHL's  from  competing  methods  are  typically  very  small. 
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-  The  consequences  of  changing  1  -  a.  The  CIE's  q»ialitatively  perform  the  same 
as  we  vary  1  -  a.  Coverage  is  sensitive  to  decreases  in  1  -  a.  For  example, 
consider  CIEj^  for  the  AR(l)  with  (p  =  0.9,  b  =  2,  and  m  =  16;  for  1  -  a  =  0.99, 
0.90,  and  0.80,  Example  5  gives  coverages  of  0.985,  0.852,  and  0.713,  resp. 

3.3.3  Summary  of  Part  I 

For  small  m,  improper  CVG's  were  usually  the  rule.  For  small  m  and  b, 

CIEj^  has  better  CVG's  than  the  other  CIE's,  while  for  small  m  and  large  b,  both 
CIEj^  and  CIE^  fared  the  best.  For  fixed  m,  high  CVG  was  most  often  accompanied 
by  high  EHL  and  SHL.  The  CIE's  performed  as  expected  by  asymptotic  theory  when 
m  and/or  b  became  large.  For  large  m  and  small  b,  all  achieved  the  desired 
coverage,  and  the  EHL's  and  SHL's  tended  to  decrease  with  increasing  d.o.f.,  as 
per  Subsection  2.3.  For  large  ro  and  b,  the  ratios  from  (14)  took  effect. 

3.4  Results  from  Part  II  of  the  Uonte  Carlo  Study 

fe  considered  the  MA(1)  (9  =  -0.9),  AR(l)  (<p  =  0.9),  EAR(l)  (tp  =  0.9),  and 
M/M/1  (p  =  0.8),  with  1  -  «  =  0.90  and  d.o.f.  d  =  3  and  15.  Table  3  gives  CVG's 
as  a  function  of  the  sample  size  n.  For  d  =  3,  CIEj^,  CIE^,  and  CIE^  perform 
about  the  same  in  terms  of  CVG;  CIE^  fares  poorly  for  small  n.  For  d  =  15  and 
small  n,  CIE^  does  a  bit  better  than  the  others  in  terms  of  CVG;  CIE^  is  not 
competitive.  However,  the  performance  of  CIE^  is  "more  variable"  than  the 
other  CIE's  over  the  range  of  stochastic  processes,  d,  and  n  under  study.  For 
instance,  the  CVG  of  CIEq  sometimes  significantly  overshoots  1  -  a,  especially 
for  small  d  (though  this  is  understandable  since  OBM  was  designed  for 
large  b) .  As  n  grows  (for  d  =  3  or  15),  it  appears  that  CIEj^,  CIE^,  and  CIE^ 
achieve  CVG  =  0.90  at  about  the  same  n. 
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4.  DISCUSSION 

We  first  discuss  the  causes  of  improper  CIE  coverage;  why  are  some  CIE's 
better  than  others?  We  then  consider  the  question  of  which  CIE  is  "best?" 

4.1  Causes  of  Improper  Coverage 

A  CIE  of  the  form  in  (1)  will  attain  perfect  coverage  if  its  associated 
pivot  T  =  (S  -  p)/(V/n)'^  ~  t(d);  this  requires  (i)  S  ~  Nor(p,a^/n),  (ii)  V  ~ 
a^X  (d)/d,  and  (iii)  independence  of  and  V.  Requirement  (i)  is  satisfied  if 
the  marginal  distribution  of  the  XVs  is  normal.  If  the  X^'s  are  not  symmetric, 
then  T  may  be  skewed  for  small  n.  But  in  most  cases,  a  central  limit  theorem 
asserts  that  (i)  approximately  holds  as  n  grows.  We  believe  that  (ii)  is  the 
key  requirement.  At  a  minimum,  V  must  be  approximately  unbiased  as  an 
estimator  of  a  (or  o^) .  In  fact,  since  variance  directly  affects  the  CIE's 
length,  we  claim  that  Bias(V)  is  often  the  main  cause  of  improper  coverage  (at 
least  for  small  m) ;  see  below.  Concerning  (iii),  Glynn  (1982)  and  Kang  and 
Goldsman  (1990)  demonstrate  that  asymmetry  in  coverage  is  directly  related  to 
dependence  between  and  V.  However,  Kang  and  Goldsman  give  examples  which 
show  that  actual  coverage  is  not  necessarily  affected  by  such  dependence.  We  do 
not  view  dependence  between  X^  and  V  as  a  direct  cause  of  improper  coverage. 

We  first  analyze  the  effect  of  Bias(V)  on  requirement  (ii)  by  examining  V  = 

a*V/E[V]  instead  of  V;  note  that  E[V]  =  -  To  illustrate,  we  shall  use  the 

NOBM  estimator  on  the  AR(1)  and  EAR(l)  processes  with  ip  =  0.9.  For  the  AR(1) 

with  b  =  2,  Example  5  says  that  ~  E[Vj^]X^(l),  and  so  a^X^(l);  for  this 

case,  correction  for  bias  results  in  precisely  the  desired  distribution. 

Empirical  p.d.f.'s  of  Vjj/a*  and  (based  on  100000  independent  runs)  are 

plotted  in  Figure  5  for  the  AR(1)  and  EAR(l)  with  b  =  8  and  various  m.  For  the 

^2  2 

AR(1),  the  sample  p.d.f.'s  of  approach  the  X  (7)  p.d.f.  as  m  increases; 
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the  corrected  is  nearly  (but  not  quite)  X  (7)  (Figures  5(a)  and  5(b)). 

The  empirical  p.d.f.'s  for  the  EAR(l)  exhibit  similar  behavior,  except  that  the 

2 

convergence  to  the  X  (7)  is  slower  (Figures  5(c)  and  5(d)).  Thus,  for  these 
examples,  correction  for  bias  mitigates  the  violation  of  (ii).  Comparing  the 
AR(l)  results  wit  -,  the  EAR(l)'s,  we  conclude  that  the  effect  of  Bias(V)  is 
greater  than  that  of  the  marginal  distribution. 

We  will  next  show  for  the  stochastic  processes  investigated  that 
correcting  for  Bias(V)  results  in  good  CIE's.  The  notation  CCVG  in  Table  1  is 
the  exact  (from  Example  5)  or  sample  coverage  obtained  for  the  AR(l)  and  EAR(l) 
models  with  ip  =  0.9  and  b  =  2  and  16.  For  the  AR(l)  with  b  =  2  and  any  m,  the 
corrected  NOBM  pivot  (1^  -  p)/(Vj^/n)‘^  ~  t(l).  So  the  CCVG's  for  ClEj^  are 
exactly  1  -  a  =  0.90;  thus,  for  this  example,  Bias(Vj^)  is  the  sole  cause  of 
improper  coverage.  The  corresponding  MC  CCVG  improvements  for  the  EAR(l)  with 
b  =  2  are  significant  but  not  as  good  as  those  for  the  AR(l),  this  indicating 
that  a  secondary  marginal  distributional  effect  is  present.  This  conclusion  is 
illustrated  yet  again  by  Figure  6.  (The  EAR(l)’s  empirical  p.d.f.'s  are 
somewhat  skewed  for  small  m  as  explained  earlier.)  The  results  from  Table  1 
and  Figures  5  and  6  suggest  that,  for  small  m,  there  are  also  tertiary 
contributors  to  improper  coverage,  perhaps  inter-batch  correlation. 

4.2  Which  is  the  best  CIE? 

The  question  of  which  CIE  is  the  "best"  depends  on  the  criteria  being 
used.  A  CIE  is  first  judged  by  the  validity  of  its  coverage.  But  the  MC  work 
showed  that  a  CIE  with  good  coverage  might  produce  relatively  wide  half- 
lengths;  so  the  E[H]  and  Var(H)  measures  can  not  be  ignored.  Since  coverage, 
E[H],  and  Var(H)  are  functions  of  the  stochastic  process,  the  CIE  in  use,  the 
number  of  batches  b,  the  batch  size  m,  and  the  level  a,  we  can  see  that  the 
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determination  of  the  "best"  CIE  is  not  straightforward. 

Results  (12)  and  (13)  say  that  E[H]  and  Var(H)  decrease  in  the  d.o.f.  d  as 
m  becomes  large  with  fixed  b;  so  the  more  d.o.f.  the  better  (although  Schmeiser 
1982  finds  that  there  is  little  to  be  gained  by  taking  d  >  30).  Nevertheless, 
it  would  be  incorrect  to  conclude  that  CIE^  (which  has  the  largest  d.o.f.  of 
those  CIE '8  under  study)  is  always  the  best.  The  Part  I  MC  results  showed  that 
for  small  b,  CIEj^  required  the  smallest  value  of  m  to  achieve  valid  coverage 
(for  large  b,  CIEj^  and  CIEq  required  the  smallest  m);  if  coverage  were  the  only 
criterion  for  CIE  comparison,  CIEj^  would  be  declared  the  best  -  not  CIE^,.  This 
shows  that  large  sample  superiority  does  not  necessarily  extend  to  the  small 
sample  case.  Indeed,  we  did  not  include  the  STS  "NOBU+maximum"  CIE  from 
Schruben  (1983),  which  has  4b-l  d.o.f.  (and  hence  superior  asymptotic 
properties),  since  it  exhibited  poor  small  sample  performance  compared  to  the 
other  CIE's  (including  CIE^).  For  instance,  for  the  AR(l)  with  ip  =  0.9,  1  -  a  = 
0.90,  b  =  2,  and  m  =  16,  64,  and  256,  the  NOBM+max  CIE  attained  CVG’s  (based  on 
1000  independent  simulation  runs)  of  0.397,  0.606,  and  0.776,  resp.;  CIE^ 
achieved  CVG's  of  0.684,  0.849,  and  0.895,  resp.  (Table  1). 

Another  basis  for  comparison  among  CIE's  is  to  determine  which  requires 
the  smallest  sample  size  n  to  achieve  valid  coverage  for  some  fixed  d.o.f.  d. 
This  was  the  aim  of  Part  II  of  the  MC  study,  where  CIEj^,  CIE^,  and  CIE^  fared 
more  or  less  the  same;  CIE^  was  not  competitive.  So  there  was  no  clear  winner 
using  the  criterion  of  coverage  under  fixed  d.o.f. 

We  can  still  make  some  recommendations  (for  fixed  sample  size 
procedures).  All  of  the  CIE's  studied  here  are  easy  to  use.  Batch  means  is 
the  simplest  method  to  understand.  All  are  asymptotically  valid  as  the  batch 
size  m  grows;  but  it  "never  hurts  too  much"  to  use  CIEj^  (in  comparison  to  the 
other  CIE's)  in  case  m  is  not  "large  enough."  The  price  to  be  paid  when  m  is 
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small  is  that  E[Hj^]  and  Var(Hj^)  are  larger  than  their  competitors,  particularly 

when  b  is  also  small.  If  the  user  is  somehow  confident  that  the  batch  sizes 

are  large  enough  to  achieve  valid  coverage,  then  the  user  should  fix  the  d.o.f. 

(perhaps  between  15  and  30  for  the  "usual"  a  values)  and  select  from  among 

CIE„,  CIE»,  and  CIE„.  However,  one  of  the  most  difficult  open  questions  in 
W  U  L 

simulation  output  analysis  is  the  determination  of  "sufficient  batch  size"  (cf. 
Fishman  1978  and  Schmeiser  and  Song  1989). 

5.  SUMMARY  AND  CONCLUSIONS 

We  studied  the  behavior  of  different  CIE's  with  special  emphasis  placed  on 
small  sample  size  performance  for  various  stationary  stochastic  processes.  If 
the  achieved  coverages  are  at  the  desired  levels,  and  if  the  total  number  of 
observations  n  is  fixed,  the  ranking  of  the  CIE’s  with  respect  to  E[H]  and 
Var(H)  is  determined  by  the  d.o.f.  each  CIE  has.  In  this  case,  the  CIE  with 
the  largest  d.o.f.  has  the  smallest  mean  and  most  stable  confidence  interval 
length. 

Perhaps  our  most  important  finding  was  that,  in  small  sample  settings,  a 
CIE  with  more  d.o.f.  may  not  actually  be  "better"  than  a  competing  CIE;  some 
CIE's  may  require  more  observations  than  others  before  the  asymptotics 
necessary  for  CIE  validity  manifest  themselves.  Quite  often,  the  CIE's  with 
the  highest  d.o.f. 's  performed  the  most  poorly  in  terms  of  coverage! 

The  bias  of  V  as  an  estimator  of  a  (or  a^)  plays  a  significant  role  in 
CIE  performance  -  the  less  bias  the  better.  For  instance,  when  m  and  b  are 
fixed,  the  relative  performance  of  the  CIE's  with  respect  to  coverage  is 
directly  related  to  Bias(V).  A  secondary  factor  in  CIE  performance  concerns 
the  underlying  marginal  distribution  of  the  X.’s.  Further,  with  fixed  m  (and 
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even  more  so  with  fixed  n) ,  coverage  often  deteriorates  as  b  or  a  increase; 

this  is  partially  due  to  the  fact  that,  in  these  cases,  t^  l-ot/2  decreases. 

Which  CIE  should  one  use  in  practice?  If  the  sample  size  n  is  "large 

enough,"  we  could  probably  argue  successfully  for  CIEj^,  CIEq,  or  CIE^  with 

common  d.o.f.,  15  <  d  <  30.  In  comparison  to  the  other  methods,  CIEj^  is 

probably  the  "safest"  small  sample  method.  There  are  several  interesting 

research  lines.  We  would  like  to  see  more  emphasis  on  small  sample  results, 

including  sequential  procedures  (which  were  not  investigated  here).  Another 

question  concerns  the  fact  that  for  fixed  d.o.f.,  CIE„,  CIE»,  and  C1E„ 

N  U  L 

achieve  approximately  valid  coverage  for  about  the  same  sample  size. 

Further,  a  good  batch  size  estimation  procedure  would  be  of  tremendous 
import . 

ACKNOWLEDGMENTS 

We  thank  M.J.  Rao  and  J.  Chacko  for  some  computer  support.  We  are  indebted 
to  Bruce  Schmeiser  for  many  interesting  discussions.  This  work  was  supported 
in  part  by  the  Air  Force  Systems  Command,  Rome  Air  Development  Center, 
Griffiss  Air  Force  Base,  New  York,  and  by  Naval  Postgraduate  School, 

Monterey,  California. 


REFERENCES 

Bratley,  P.,  B.L.  Fox,  and  L.E.  Schrage.  1987.  A  Guide  to  Simulation , 

2nd  Ed.,  Spr inger-Ver lag.  New  York. 

Conway,  R.W.  1963.  "Some  Tactical  Problems  in  Digital  Simulation." 
Management  Science,  10,  47-61. 

Crane,  M.A.  and  D.L.  Iglehart.  1975.  "Simulating  Stable  Stochastic  Systems: 
III.  Regenerative  Processes  and  Discrete-Event  Simulations."  Operations 
Research,  23,  33-45. 

Crane,  M.A.  and  A.J.  Lemoine.  1977.  An  Introduction  to  the  Regenerative 
Method  for  Simulation  Analysis.  Spr inger-Ver lag,  Berlin. 

Damerdji,  H.  1987.  Topics  in  Discrete-Event  Stochast ic  Systems.  Ph.D. 
Dissertation,  Dept,  of  Industrial  Engineering,  Univ.  of  Wisconsin,  Madison, 
Wis. 

Fishman,  G.S.  1973.  Concepts  and  Methods  in  Discrete  Event  Digital 
Simulation.  John  Wiley  A  Sons,  New  York. 


-18- 


Fishman,  G.S.  1978.  Principles  of  Discrete  Event  Simulation.  Wiley,  NY. 

Glynn,  P.W.  1982.  "Coverage  Error  for  Confidence  Intervals  Arising  in 
Simulation  Output  Analysis."  Proc.  1982  Winter  Simul.  Conf.,  369-375. 

Glynn,  P.W.  and  D.L.  Iglehart.  1990.  "Simulation  Output  Analysis  Using 
Standardized  Time  Series."  Mathematics  of  Operations  Research,  15,  1-16. 

Goldsman,  D.  1984.  On  Using  Standardized  Time  Series  to  Analyze  Stochastic 
Processes.  Ph.D.  Dissertation,  School  of  Operations  Research  and  Industrial 
Engineering,  Cornell  Univ.,  Ithaca,  NY. 

Goldsman,  D.,  K.  Kang,  and  R.G.  Sargent.  1986.  "Large  and  Small  Sample 
Comparisons  of  Various  Variance  Estimators."  Proc.  1986  Winter  Simul.  Conf., 
278-284. 

Goldsman,  D.  and  M.S.  Ueketon.  1986.  "A  Comparison  of  Several  Variance 
Estimators."  Tech.  Report  jjfJ-85-12,  School  of  Industrial  and  Systems 
Engineering,  Georgia  Institute  of  Technology,  Atlanta. 

Goldsman,  D.  and  L.  Schruben.  1984.  "Asymptotic  Properties  of  Some  Confidence 
Interval  Estimators  for  Simulation  Output."  Management  Science,  30,  1217-1225. 

Heidelberger ,  P.  and  P.D.  Welch.  1981.  "A  Spectral  Method  for  Confidence 
Interval  Generation  and  Run  Length  Control  in  Simulations."  Comm.  ACM,  24, 
233-245. 

Heidelberger,  P.  and  P.D.  Welch.  1983.  "Simulation  Run  Length  Control  in  the 
Presence  of  an  Initial  Transient."  Operations  Research,  31,  1109-1144. 

Kang,  K.  1984.  Confidence  Interval  Estimation  via  Batch  Means  and  Time  Series 
Modeling.  Ph.D.  Dissertation,  School  of  Industrial  Engineering,  Purdue  Univ., 

W.  Lafayette,  IN. 

Kang,  K.  and  D.  Goldsman.  1990.  "The  Correlation  between  Mean  and  Variance 
Estimators  in  Computer  Simulations."  HE  Transactions,  22,  15-23. 

Law,  A.M.  and  W.D.  Kelton.  1984.  "Confidence  Intervals  for  Steady-State 
Simulations,  I:  A  Survey  of  Fixed  Sample  Size  Procedures."  Operations 
Research ,  32,  1221—1239. 

Lewis,  P.A.W.  1980.  "Simple  Models  for  Positive-Valued  and  Discrete-Valued  Time 
Series  with  ARMA  Correlation  Structure."  Multivariate  Analysis-V, 

P.R.  Krishnaiah,  ed..  North  Holland  Publishing  Company. 

Meketon,  M.S.  1980.  The  Variance  Time  Curve:  Theory,  Estimation  and 
Application.  Ph.D.  Dissertation,  School  of  Operations  Research  and  Industrial 
Engineering,  Cornell  Univ.,  Ithaca,  NY. 


Meketon,  M.S.  and  B.W.  Schmeiser.  1984.  "Overlapping  Batch  Means:  Something 
for  Nothing?"  Proc.  1984  Winter  Simul.  Conf.,  227-230. 


-19- 


Uoran,  P.A.P.  1975.  "The  Estimation  of  Standard  Errors  in  Monte  Carlo 
Simulation  Experiments."  Biometrika,  62,  1-4. 

Muirhead,  R.J.  1982.  Aspects  of  Multivariate  Statistical  Theory.  Wiley,  NY. 

Priestley,  M.B.  1982.  Spectral  Analysis  and  Time  Series.  Academic  Press,  NY. 

Sargent,  R.G.,  K.  Kang,  and  D.  Goldsman.  1987.  "An  Investigation  of  Small 
Sample  Behavior  of  Confidence  Interval  Estimation  Procedures."  Working  Paper 
#87-005,  Dept,  of  Industrial  Engineering  and  Operations  Research,  Syracuse 
Univ.,  Syracuse,  NY  13244. 

Schmeiser,  B.W.  1982.  "Batch  Size  Effects  in  the  Analysis  of  Simulation 
Output."  Operations  Research,  30,  556-568. 

Schmeiser,  B.W.  1986.  Personal  Communication. 

Schmeiser,  B.W.  and  W.T.  Song.  1989.  "Optimal  Mean-Squared-Error  Batch  Sizes." 
Tech.  Report,  School  of  Industrial  Engineering,  Purdue  Univ.,  W.  Lafayette,  IN. 

Schriber,  T.J.  and  R.W.  Andrews.  1984.  "ARMA-Based  Confidence  Intervals  for 
Simulation  Output  Analysis."  Amer.  J.  Math,  and  Mgmt.  Sci.,  4,  345-373. 

Schruben,  L.  1983.  "Confidence  Interval  Estimation  Using  Standardized  Time 
Series."  Operations  Research,  31,  1090-1108. 

Song,  W.T.  and  B.W.  Schmeiser.  1989.  "Estimating  Variance  of  the  Sample  Mean: 
Quadratic  Forms  and  Cross-Product  Covariances."  Operations  Research  Letters, 

7,  259-266. 


APPENDIX 


Proof  of  (2)-(5):  We  will  use  the  following  facts: 


j;k  i  .  5;l<  ipl  =  p[l-(;ktl^pk4kpk->l: 

=  l  1-p  il=l 

For  a  covariance  stationary  process  with  =  Cov(X^ ,X^^j^) , 

=  n-Var(X  )  =  Y  +  -  Y”  |  (n-i)Y.. 
n  n'^  0  n^i=l  i 


(A-1) 


(A-2) 
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Since 

.  Z5-l  -  •’iE[S^-B[5j]i  =  b[Var(X„)-Var(S_,)].  (A-i) 

result  (2)  follows  by  (A-3)  and  simplification  tsee  Moran  1975). 


We 


also  have  E[Vq]  =  E-lf  ^  E[(X(i.m)  -  1^)^] 


=  (n-m+l)[Var(X^)  +  Var(X^)]  -  E[a(  i  .m)Xj^] . 


(A-5) 


Since  p  =  0.  mnE[X(i.m)X^]  =  j;““J  E[X.^.Xj^]  =  "  Xk=l  ^j-k 


m+i-l  T>n 


;am+i-l  J-k  X  vi'  k-j^  ^  r.m+i-1  1  t  y  -  >f\- 

k=j  +  l  ^  ^  ‘^j^i  l-ip 


=  Zjir'  ai-.,  ^  z 


=  m 


l+(p  ma/  i  ,  n-m+2  -ia  a2 

’T^  -  (1-<P  )^‘P  +  ‘P  “P  )/(i-<p)  • 


(A-7) 


We  obtain  (3)  by  substituting  (A-3)  and  (A-7)  into  (A-5).  , 


// 


(4)  follows  from  (A. 5-15)  of  Goldsman  (1984);  (5)  follows  from  (2)  and  (4). 


// 


Proof  of  (6)-(9):  The  MA(1)  has  covariance  function  =  1  +  9^,  =  9,  and 

Yj_  =  0,  otherwise.  By  (A-2), 

-  (1+9)^  -  29/n  and  =  lim^_^^0^  =  (1+9)^.  (A-8) 

Result  (6)  then  follows  from  (A-4)  and  (A-8). 


Note  that 


yn 

^k=l 


j-k 


■  (1+9)^ 

.  (l+9)^-9 


if  1  <  j  <  n 


if  j  =  1  or  n 
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So  by  (A-6) , 

m(l+©l^  if  1  <  i  <  n-m+1 

mnE[X(i,m)X  ]  =  '  2  ’  (A-9) 

in(l+6)  -0  if  j  =  1  or  n-m+1 


We  obtain  (7)  by  substituting  (A-8)  and  (A-9)  into  (A-5). 


(8)  follows  from  (A. 5-23)  of  Goldsman  (1984);  (S)  follows  from  (6)  and  (8). 


. ,X^),  where  the 


X^'s  are  stationary  normal  with  covariance  matrix  Z.  Suppose  G  is  an  n x n 
symmetric  matrix  and  V  =  X'GX.  Problem  1.22  of  Muirhead  (1982)  says  that  V  and 


X^  are  independent  iff  I'SG  =  O',  where  1'  [O']  is  a  Ixn  vector  of  1  s  [O's]. 

Song  and  Schmeiser  (1989)  note  that  V„  =  X'G^X  and  V  =  X'G.X,  where 

N  ~  N~  A  A~ 


J  -J 
-J  J 


=  (hh.)  for  b  =  1, 
1  J 


where  J  is  an  n/2 x n/2  matrix  consisting  o+  I's,  and  h^  =  (n+l)/2  -  i, 

1  =  1,  ..,n.  Then  Gj^  and  G^  both  meet  Muirhead's  condition  for  any  E.  i. 
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Table  1  -  Some  Results  for  the  AR(l)  and  EAR(l)  Uodels  with  (p  =  0.9.  All 
entries  for  E[V]  are  exact.  Results  for  CVG  and  CCVG  arOgexactC)  or  are  based 
on  2000  independent  simulation  runs.  For  these  models,  o  =  19. 

E[V]  AR(1)  EAR(l) 


(or  o^)  b  =  2  b  =  16  b  =  2  b  =  16 

^  n 


b=2 

b=16 

CVG 

CCVG 

CVG 

CCVG 

CVG 

CCVG 

CVG 

CCVG 

m  =  4 

NOBU 

0.86 

2.68 

0.744* 

0.900* 

0.479 

0.866 

0.625 

0.798 

0.465 

0.810 

OBU 

0.58 

2.64 

0.453 

0.840 

0.466 

0.859 

0.337 

0.673 

0.446 

0.804 

Area 

0.31 

0.31 

0.418 

0.907 

0.183 

0.888 

0.296 

0.701 

0.164 

0.799 

Comb 

0.49 

1.46 

0.432 

0.878 

0.365 

0.632 

0.321 

0.683 

0.337 

0.806 

^  2 
^  n 

6.19 

16.19) 

m  =  16 

NOBM 

6. 10 

9.27 

0.852* 

0.900* 

0.756 

0.896 

0.791 

0.850 

0.742 

0.872 

OBM 

4.22 

9.23 

0.673 

0.864 

0.750 

0.890 

0.609 

0.^73 

0.733 

0.866 

Area 

2.85 

2.85 

0.684 

0.896 

0.491 

0.891 

0.615 

0.800 

0.480 

0.866 

Comb 

3.94 

5.96 

U.689 

0.884 

0.651 

0.893 

0.625 

0.795 

0.646 

0.868 

/  2 
(<7 
^  n 

13.57 

18.30) 

m  =  64 

NOBM 

14.79 

16.02 

0.891* 

0.900* 

0.874 

0.904 

0.881 

0.389 

0.869 

0.892 

OBM 

12.84 

16.00 

0.841 

0.883 

0.862 

0.894 

0.810 

0.841 

0.867 

0.888 

Area 

11.29 

11.29 

0.849 

0.897 

0.801 

0.894 

0.822 

0.866 

0.803 

0.884 

Comb 

12.45 

13.58 

0.852 

0.895 

0.840 

0.898 

0.819 

0.858 

0.841 

0.888 

(  2 

17.59 

18.82) 

ro  =  256 

NOBM 

17.95 

18.25 

0.898* 

0.900* 

0.899 

0.905 

0.900 

0.902 

0.886 

0.890 

OBM 

17.30 

18.25 

0.895 

0.904 

0.902 

0.905 

0.900 

0.909 

0.890 

0.897 

Area 

16.90 

16.90 

0.895 

0.902 

0.886 

0.903 

0.888 

0.895 

0.878 

0.896 

Comb 

17.25 

17.56 

0.890 

0.901 

0.897 

0.908 

0.896 

0.902 

0.880 

0.892 

^  2 

18.65 

18.96) 

m  =  1024 

NOBM 

18.74 

18.81 

0.900* 

0.900* 

0.897 

0.899 

0.883 

0.883 

0.909 

0.910 

OBM 

18.56 

18.81 

0.905 

0.907 

0.897 

0.898 

0.915 

0.918 

0.902 

0.903 

Area 

18.47 

18.47 

0.901 

0.905 

0.893 

0.897 

0.893 

0.895 

0.898 

0.901 

Comb 

18.56 

18.64 

0.892 

0.896 

0.900 

0.902 

0.895 

0.897 

0.899 

0.903 

(a^  18.91  18.99) 

n 
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Table  2  -  CVG  Results  for  the  AR(1)  Process,  b  =  16,  1  -  a  =  0.90. 
All  entries  are  based  on  at  least  1000  independent  simulation  runs. 


<p 

-0.9 

-0.5 

m  =  4 

NOBW 

0.962 

0.945 

OBU 

0.966 

0.949 

Area 

1.000 

0.974 

Comb 

0.998 

0.971 

m  =  16 

NOBM 

0.951 

0.915 

OBU 

0.958 

0.922 

Area 

0.988 

0.942 

Comb 

0.981 

0.926 

m  =  64 

NOBM 

0.918 

0.904 

OBU 

0.929 

0.916 

Area 

0.943 

0.901 

Comb 

0.931 

0.900 

m  =  256 

NOBM 

0.915 

0.903 

OBM 

0.915 

0.911 

Area 

0.931 

0.914 

Comb 

0.921 

0.918 

m  =  1024 

NOBM 

0.906 

0.901 

OBM 

0.903 

0.899 

Area 

0.907 

0.904 

Comb 

0.909 

0.906 

0.0  0.5  0.9 


0.913  0.823  0.479 
0.908  0.816  0.466 
0.898  0.677  0.183 
0.900  0.755  0.365 


0.905  0.887  0.756 
0.906  0.895  0.750 
0.908  0.860  0.491 
0.910  0.880  0.651 


0.903  0.901  0.874 
0.910  0.908  0.862 
0.885  0.885  0.801 
0.896  0.894  0.840 


0.904  0.900  0.899 
0.909  0.909  0.902 
0.918  0.917  0.886 
0.916  0.915  0.897 


0.899  0.899  0.897 
0.897  0.897  0.897 
0.906  0.904  0.893 
0.902  0.900  0.900 
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Table  3  -  CVG  Results  for  Variance  Estimators  with  Common  d.o.f. 
All  entries  are  based  on  2000  independent  simulation  runs. 


n  = 

32 

64 

128 

256 

512 

1024 

2048 

4096 

8192 

AR(1). 

NOBM 

<0  =  0.9 
0.688 

.  d  =  3 
0.787 

0.848 

0.880 

0.889 

0.899 

0.896 

0.902 

0.905 

OBU 

0.692 

0.792 

0.834 

0.874 

0.896 

0.910 

0.903 

0.923 

0.920 

Area 

0.505 

0.674 

0.795 

0.858 

0.880 

0.900 

0.895 

0.889 

0.912 

0.701 

0.799 

0.847 

0.886 

0.883 

0.903 

0.895 

0.905 

0.908 

EARfl) 

NOBM 

.<0  =  0- 
0.622 

9.  d  = 
0.736 

3 

0.826 

0.862 

0.887 

0.900 

0.896 

0.903 

0.896 

OBM 

0.631 

0.725 

0.820 

0.864 

0.887 

0.906 

0.917 

0.912 

0.916 

Area 

0.426 

0.630 

0.771 

0.844 

0.868 

0.886 

0.893 

0.893 

0.893 

0.636 

0.740 

0.822 

0.863 

0.881 

0.898 

0.894 

0.902 

0.896 

M/M/1, 

NOBM 

D  =  0.8 
0.437 

.  d  =  3 
0.541 

0.632 

0.707 

0.768 

0.804 

0.836 

0.862 

0.877 

OBM 

0  432 

0.540 

0.627 

0.704 

0.752 

0.799 

0.842 

0.867 

0.886 

Area 

0.290 

0.420 

0.533 

0.624 

0.708 

0.769 

0.817 

0.835 

0.855 

0.442 

0.548 

0.634 

0.701 

0.758 

0.806 

0.838 

0.860 

0.871 

MA.(l), 

NOBM 

e  =  -0. 
0.985 

9.  d  = 
0.981 

3 

0.986 

0.975 

0.953 

0.932 

0.931 

0.920 

0.911 

OBM 

1.000 

1.000 

1.000 

1.000 

0.998 

0.986 

0.971 

0.959 

0.944 

Area 

0.994 

0.993 

0.989 

0.980 

0.967 

0.955 

0.947 

0.933 

0.918 

Comb 

0.994 
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Figure  2:  AR(l)  Process  Performance  Measures 
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Figure  4:  M/M/1  Process  Performance  Measures 
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Figure  6:  Empirical  p.d.f.'s  (based  on  100000  runs)  of  (X^  -  n)/(Vj^/n)^  and 
(X^  -  for  the  AR(l)  and  EAR(l)  Kith  (p  =  0.9  and  b  =  0 


Distribution  List 


Defense  Technical  Information  Center 
Cameron  Station 
Alexandria,  VA  22314 

Dudley  Knox  Library,  Code  0142 
Naval  Postgraduate  School 
Monterey,  CA  93943 

Office  of  Research  Administration 
Code  08 

Naval  Postgraduate  School 
Monterey,  CA  93943 

Library,  Center  for  Naval  Analyses 
4401  Ford  Avenue 
Alexandria,  VA  22302-0268 

Department  of  Administrative  Sciences  Library 
Code  AS 

Naval  Postgraduate  School 
Monterey,  CA  93943 

Prof.  Keebom  Kang 
Code  AS/Kk 

Administrative  Sciences  Department 
Naval  Postgraduate  School 
Monterey,  Ca.  93943 

Prof.  Robert  G.  Sargent 
Simulation  Research  Group 
Syracuse  University 
Syracuse,  NY  13244 

Prof.  David  Goldsman 
School  of  Industrial  &  Systems  Engineering 
Georgia  Institute  of  Technology 
Atlanta,  GA  30332 


1 


1 


10 


20 


10 


